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Chapter 15. Modeling of Data 



15.0 Introduction 



Given a set of observations, one often wants to condense and summarize the 
data by fitting it to a "model" that depends on adjustable parameters. Sometimes the 
model is simply a convenient class of functions, such as polynomials or Gaussians, 
and the fit supplies the appropriate coefficients. Other times, the model's parameters 
come from some underlying theory that the data are supposed to satisfy; examples 
are coefficients of rate equations in a complex network of chemical reactions, or 
orbital elements of a binary star. Modeling can also be used as a kind of constrained, 
interpolation, where you want to extend a few data points into a continuous function, 
but with some underlying idea of what that function should look like. 

The basic approach in all cases is usually the same: You choose or design a 
figure-of-merit Junction C*merit function,** for short) that measures the agreement 
between the data and the model with a particular choice of parameters. The merit 
function is conventionally arranged so mat small values represent close agreement. 
The parameters of the model are then adjusted to achieve a rnimmum in the merit 
function, yielding best-fit parameters. The adjustment process is thus a problem in 
miriimization in many dimensions. This optimization was the subject of Chapter 10; 
however, there exist special, more efficient, methods that are specific to modeling, 
and we will discuss these in mis chapter. 

There are important issues that go beyond the mere finding of best-fit parameters. 
Data are generally not exact. They are subject to measurement errors (called noise 
in the context of signal-processing). Thus, typical data never exactly fit the model 
that is being used, even when that model is correct We need the means to assess 
whether or not the model is appropriate, that is, we need to test the goodness-of-fit 
against some useful statistical standard. 

We usually also need to know the accuracy with which parameters are de- 
termined by the data set. In other words, we need to know the likely errors of 
the best-fit parameters. 

Finally, it is not uncommon in fitting data to discoverjhat the merit function 
is not unimodal, with a single rninimum. In some cases, we may be interested in 
global rather than local questions. Not, "how goodTs this fit?" but rather, "how sure 
am I that there is not a very much better fit in some comer of parameter space?" 
As we have seen in Chapter 10, especially §10.9, this kind of problem is generally 
quite difficult to solve. 

The important message we want to deliver is that fitting of parameters is not 
the end-all of parameter estimation. To be genuinely useful, a fitting procedure 
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should provide (i) parameters, GO error estimates on the parameters, and (Hi) a 
sta t i stic al measure of goodness-of-fit. When the third item suggests that the model 
is an unlikely match to the data, then items (i) and (n) are probably worthless. 
Unfortunately, many' practitioners of parameter estimation never proceed beyond 
item (i). They deem a fit acceptable if a graph of data and model "looks good.'* This 
approach is known as chi-by-eye. Luckily, its practitioners get what they deserve. 

crre> references and further reading: 

Bevington, PSL 1969, Data Reduction and Error Analysis for tte Physical Sciences (New York: 
McGraw-Hnf). 

Brovmleo, KK 1965. Statistical Theory and Methodology t 2nd ect (New York: YVBey). 
Martin. BA 1 971, Statistics for PhysJ&sts (New York: Academic Press), 
von Msgs, R. 1964, Mathematical Theory of ProbabtSty and Statistics (New Yortc Academic 
Press), Chapter X. 

Konv G-A^ and Kom, T.ML 1968, Matherna^cairUaidbookfofSda^tsaiidEn^neerSt 2nd ed. 
(New York: McGraw-HO), Chapters 18-19. 



15.1 Least Squares as a Maximum Likelihood 
Estimator 

Suppose that we are fitting N data points (x$ t t = 7V\ to a nrodel that 
has M adjustable parameters Oj 9 j ~ 1,... ,Af. The model predicts a functional 
relationship between the measured independent and dependent variables, 

= y(x; ai . . , a M ) (15.1.1) 

where the dependence on the pararneters is indicated explicitly on the right-hand side. 

What, exactly, do we want to mmiim» to get fitted values for the Oj *s? The 
first thing that comes to mind is the famiHar least-squares fit, 

H 

rhhnxnize overai . . . aj/ ; ^ fefe — Qq . . ■ ajflf (15.1.2) 

But where does this come from? What general principles is it based on? The answer 
to these questions takes us into the subject of maximum likelihood estimators. 

Given a particular data set of x»*s and k's, we have the intuitive feeling that 
some parameter sets a* ... a are very unHkery — those for which the model 
function y(x) looks nothing like the data — while others may be very likely — those 
that closely resemble the data. How can we quantify this intuitive reeling? How can 
we select fitted parameters that are "most likely" to be correct? It is not nieamngful 
to ask the question, "What is the probability that a particular set of fitted parameters 
«i - - • <*>m is correct?** The reason is that there is no statistical universe of models 
from which the parameters are drawn. There is just one model* the correct one, and 
a statistical universe of data sets that are drawn from ill 
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■mat beine the case, we can, however, turn the question around, and ask, "Given 
a \7f parsers, what is the probability that this data setcouWfcwe 

If Ae y^s W on continuous values, the probability will ^ ^ 
0< * unea ' " .t! ^ " olus or mirms some fixed Ay on each data point. 
tZZ:TyZZ^^™ ^rstood, ff^probabrntyofobt^the 
Stasetts^itesimali/small. then we can conclude that the patten »unde* 
Edition axe unlikely" to be right. Conversely, our ** °* 

data set should not be too improbable for the correct choice ofP^^" 

In other words, we identify the probability of the data pven the parameters 
fwhiS is alTh^maticaUy computable number), as the livelihood of the parameters 
SSfto daT-rhis identification is entirely based on It ha, rK, Tc^al 

mXmatical basis in and of itself; as we already remarked, statics is not a 

^tcl intuiave identification, however His jjj. -jJWbj 

««, to decide to fit for the parameters 0l . . . <*m precisely by finding those values 
2 L liSihocTo^ed in the above way. This form of parameter 

estimation is maximum likelihood estimation. c™™^ that each data 

w . __ now ready to make the connection to (15.1-2). Suppose tnat eacn oaia 
point^£ IZ^Lv* error that is independently ^ - d 
nZ^clsian) distribution around the "true" model y(*). 
^d^j deviations a of these normal distributions are the same for aQ points lien 
r^bS offte data set is tne product of ti* probabilities of each point. 

Notice that mere is a factor Ay in each term in me product. Marirnizing j^"*^** 
Lqu^dcnt to maxirmzing its logarithm, or minirnizing the negatrve of its loganthm. 

namely. 

Since N, c, and Ay are a* constants, mirnmizing this equation is equivalent to 

"'"^Ji Aat least-squares fining is a rnaximuffi^elmoo^estirn^ 
of thHnU^meters 4 me me^ement errors are ir^ependent and r^y 
^striLX^nstant standard deviation. Notice that we made no action 
„iw™t the linearitv or nontinearity of the model vfo m • -) m its parameters 
SeToXwTwiU rel« our assumption of constant standard demons 
^d ob^n the very similar formulas for what fa caBed "enr-sq^are fimng ^ 
^rijhted least-squares fitting." first, however, let us discuss further our very 

a^TdXl 3 ZTJSS statisticians have been in love w^ 
^ faTtLthe proclbuity distribution of me sum of a very large n^bexof very 
smaU random delations almost always converges to a normal ^mbuoon. (For 
t^^nen^ of this cenrr*/ Umit theorem, consult HI or mer 
oTmatiKmatical statistics.) Tins infatuation tended to focus mterest away from the 
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fact that, for real data, the normal distribution is often rather poorly realized, if it is 
realized at alL We are often taught, rather casually, that, on average, measurements 
will fall within db<r of the true value 68 percent of the time, within ±2o 95 percent 
of the time, and within ±3a 99.7 percent of the tone. Extending this, one would 
expect a measurement to be off by ±20a only one time out of 2 x 10 88 . We all 
know that "glitches" are much more likely than that\ 

In some instances, the deviations from a normal distribution are easy to 
understand and quantify. For example, in measurements obtained by counting 
events, the measmernent errors are usually distributed as a Poisson distribution, 
whose cumulative probability function was already discussed in §6.2. When the 
number of counts going into one data point is large, the Poisson distribution converges 
towards a Gaussian, However, the convergence is not uniform when measured in 
fractional accuracy. The more standard deviations out on the tail of the distribution, 
the larger the number of counts must be before a value close to the Gaussian is 
realized. The sign of the effect is always the same: The Gaussian predicts that "tail** 
events are much less nicely than they actually (by Poisson) are. This causes such 
events, when they occur, to skew a least-squares fit much more than they ought. 

Other frww^*, the deviations from a normal distribution are not so easy to 
understand in detail. Experimental points are occasionally just way off. Perhaps 
the power flickered during a point's measurement, or someone kicked the apparatus, 
or someone wrote down a wrong number. Points like this are called outliers. 
They can easily turn a least-squares fit on otherwise adequate data into nonsense. 
Their probability of occurrence in the assumed Gaussian model is so small that the 
maximum likelihood estimator is willing to distort the whole curve to try to bring 
them, mistakenly, into line. 

The subject of robust statistics deals with cases where the normal or Gaussian 
model is a bad approximatiori, or cases where outliers are important. We will discuss 
robust methods briefly in §15.7. AH the sections between this one and that one 
assume, one way or the other, a Gaussian model for the measurement errors in the 
data. It it quite important that you keep the limitations of that model in mind, even 
as you use the very useful methods that follow from assuming it. 

Finally, note that our discussion of measurement errors has been limited to 
statistical errors, the kind that will average away if we only take enough data. 
Measurements are also susceptible to systematic errors that will not go away with 
any amount of averaging. For example, the calibration of a metal meter stick might 
depend on its temperature. If we take all our measurements at the same wrong 
temperature, then no amount of averaging or numerical processing will correct for 
this uiirecognized systematic error. 



ChfSquare Fitting 

We considered the chi-square statistic once before, in §143. Here it arises 
in a slightly different context. 

If each data point jft) has its own, known standard deviation &u then 
equation (15.13) is modified only by putting a subscript i on the symbol <r. That 
subscript also propagates docilely into (15.1.4), so that the maximum likelihood 
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estimate of the model parameters is obtained by minimizing the quantity 

called the "chi-square." 

To whatever extent the measurement errors actually are normally distributed, 
the quantity x* is correspondingly a sum of JV squares of normally distributed 
quantities, each normalized to unit v ariance. Once we have adjusted the a x ... a M to 
minimize the value of x 2 . the terms in the sum are not all statistically independent. 
For models that are linear in the a's, however, it turns out that the probability 
distribution for different values of x 2 at its nimimum can nevertheless be derived 
analytically, and is the chi-square distribution forN-M degrees of freedom. 
We learned how to compute this probability function using the incomplete gamma 
function gaimnq in §6.2. In particular, equation (6.2.18) gives the probability Q that 
the chi-square should exceed a particular value x 3 by chance, where v = 7V - M is 
the number of degrees of freedom. The quantity <?, or its complement P = 1 - Q, is 
frequently tabulated in appendices to statistics books, but we generally find it easier 
to use gammq and compute our own values; Q = gammq (0-5x/, 0.5x 2 )- It is quite 
common, and usually not too wrong, to assume that the chi-square distribution holds 
even for models that are not strictly linear in the a's. 

This computed probability gives a quantitative measure for the goodness-of-fit 
of the model- If Q is a very small probability for some particular data set, then the 
apparent discrepancies are unlikely to be chance fluctuations. Much more probably 
either CO the model is wrong — can be statistically rejected, or (ii) someone has lied to 
you about the size of the measurement errors ^ — they are leaUy larger than stated. 

It is an important point that the chi-square probability Q does not directly 
measure the credibility of the assumption that the measurement errors are normally 
distributed. It assumes they are. In most, but not all, cases, however, the effect of 
nonnormal errors is to create an abundance of outlier points. These decrease the 
probability Q 9 so that we can add another possible, though less defmitive, conclusion 
to the above list: Qu) the measurement errors may not be normally distributed. 

Possibility (iii) is fairly conimon, and also fairly benign. It is for this reason 
that reasonable experimenters are often rather tolerant of low probabilities Q. It is 
iKHuncoinmontodeemacceptaW 0.001. 
This is not as sloppy as it sounds: Truly wrong models will often be rejected with 
vastly smaller values of Q, 10~ w , say. However, if day-in and day-out you find 
yourself accepting models with Q ~ 10~ 3 , you really should track down the cause. 

If you happen to know the actual distribution law of your measurement errors, 
then you might wish to Monte Carlo simulate some data sets drawn from a particular 
model, cf. §7^-573. You can then subject these synthetic data sets to your actual 
fitting procedure, so as to determine both the probability distribution of the x 
statistic, and also the accuracy with which your model parameters are reproduced 
by the fit- We discuss this further in § 15.6V The technique is very general, but it 
can also be very expensive. _ 

At the opposite extreme, it sometimes happens that the probab&tfy Q is too large, 
too near to 1, literally too good to be true! Nonnormal measurement errors cannot 
in general produce this disease, since the normal distribution is about as -compact" 



SDOCIO. <XP 219106BA_I_> 



15J2 Fitting Data to a Straight Une 



661 



as a distribution can be Almost always, the cause of too good a dri-square fit 
is that the experimenter, in a "fit** of conservativism, has overestimated his or her 
measurement errors. Very rarely, too good a chi-square signals actual fraud, data 
mat has been "fudged?* to fit the modeL 

A rale of thumb h mat a "typical" value of x 2 for a "moderately" good fit is 
X 2 More precise is me statement that the ^statistic 
deviation y/2*/ y and, asymptotically for large f, becomes normally distributed. 

In some cases me uncertainties associated with a set of measurements are not 
known in advance, and considerations related to x 2 fitting are used to derive a value 
for a. If we assume that all measurements have the same standard deviation, a. — o\ 
and that the model does fit well, then we can proceed by first assigning an arbitrary 
constant a to all points, next fitting for the model parameters by minimizing x 2 » 
and finally recomputing 

N 

o 2 - jy* ~ - M ) (15.1.6) 

Obviously, this approach prohibits an independent assessment of goodness-of-fit, a 
fact occasionally missed by its adherents. When, however, the measurement error 
is not known, this approach at least allows some kind of error bar to be assigned 
to the points. 

If we take the derivative of equation (15.1.5) with respect to the parameters a kt 
we obtain equations that must hold at the chi-square minimum, 

°=g(^)( £!fe ie- i ) * <»•■•" 

Equation ( 15. 1 .7) is, in general, a set of Af nonlinear equations for the Af unknown 
a*. Various of the procedures described subsequently in tins chapter derive from 
(15.1.7) and its specializations. 

CITED REFERENCES AND FURTHER REACHNG: 

Bevftngton, PJL 1969, Data Reduction and EnwAnafy& tar the Physk& Sciences Q^ybnc 

McGraw+BB). Chapters 1-4. 
von Mtees, R. 1964. Mathematical Theory of PnsbabBHy and Statistics (New \brfc Academic 

Press). §VLC. [11 



15.2 Fitting Data to a Straight Une 

A concrete example win make the considerations of the previous section more 
meaningful. We consider the problem of fitting a set of TV* data points (x f , y,-) t 
a straight-line model 



(15.2.1) 
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This problem is often called linear regression, a terminology that originated, long 
ago, in the social sciences. We assume thai the uncertainty a* associated with 
each measurement y; is known, and that the x t *s (values of the dependent variable) 
are known exactly. 

To measure how well the model agrees with the data, we use the chi-square 
merit function (15.1.5), which in this case is 

If the measurement am are normally distributed, then this merit function 
maximum likelihood parameter estimations of a and Ir. if the errors arc not ^naUy 
distributed, then the estimations are not maximum likelihood, but may soil be useful 
in a practical sense. In §15.7. we will treat the case where outlier points are so 
numerous as to lender the x 2 merit function useless. . . ,,_ 

Equation (15.2.2) is niinimized to determine a and 6. At its minimum, 
derivatives of ^(o.b) with respect to o, 6 vanish. 

da £i (15.2.3) 

These conditions can be rewritten in a convenient form if we define the following 

sums; 

tZ*i i=»°« 1=1 • (15.2.4) 

With these definitions (15.23) becomes 

The solution of these two equations in two unknowns is calculated as 
A = SS«-(S*) 2 



o = 



S„S V - 5,gr» (15.2.6) 



6 - A 



Equation (15.2.6) gives the solution for the best-fit model parameters a and 6. 
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We are not done, however. We must estimate the probable uncertainties in 
the estimates of a and 6, since obviously the measurement errors in the data must 
mtroduce some nncertamty in the determination of those parameters. If the data 
are independent, then each contributes its own bit of uncertainty to the parameters. 
Consideration of propagation of errors shows that the variance <r? in the value of 
any function will be 

For the straight line, the derivatives of a and b with respect to y t can be directly 
evaluated from the solution: 

da S„- S x Xj 
tyi of A 

db _ s Xi -s x C* 5 - 2 -*) 

Summing over the points as in (15.2.7), we get 

which are the variances in the estimates of a and 6, respectively. We win see in 
§15.6 that an additional number is also needed to characterize properry the probable 
micertainty of the parameter estimation. That number is the covariance of a and 6, 
and (as we will see below) is given by 

Coy(a, 6) = -&/A (15.2.10) 

The coefficient of correlation between the uncertainty in a and the micertainty 
in 6, which is a mtmber between -1 and 1, follows from (15.2.10) (compare 
equation 14.S.1X 

A positive value of r«* indicates that the errors in a and 6 are likely to have the 
same sign, while a negative value indicates the errors are anncorrelated, likely to 
have opposite signs. 

We are still not done. We must estimate the goodness-of-fit of the data to the 
modeL Absent this estimate, we have not the slightest mdicathm that the parameters 
a and b in the model have any meaning at all! The probability Q. that a value of 
cm-square as poor as the value (15JL2) should occur by chance is 



(15.2.12) 
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Here gammq is our routine for the incomplete gamma function Q(a, x), §6.2. If 
Q is larger than, say, 0.1, then the goodness-of-fit is believable. If it is larger 
than, say, 0.001, then the fit may be acceptable if the errors are nonnormal or have 
been moderately underestimated. If Q is less than 0.001 men the model and/or 
estimation procedure can rightly be called into question. In this latter case, turn 
to §15.7 to proceed further. 

If you do not know the individual measurement errors of the points a u and are 
proceeding (dangerously) to use equation (15.1.6) for estimating these errors, then 
here is the procedure for estimating the probable uncertainties of the parameters a 
and 6: Set a» = 1 in all equations through (15.2.6), and multiply cr 0 and a by as 
obtained from equation (15-2.9), by the additional factor y/^{{N -2), where x 
is computed by (15.2.2) using the fitted parameters a and 6. As discussed above, 
this procedure is equivalent to assuming a good fit, so you get no independent 
goodness-of-fit probability Q. 

In §143 we promised a relation between the linear correlation coefficient 
r (equation 143.1) and a goodness-of-fit measure, x 2 (equation 15^2). For 
unweighted data (all Oi — 1), that relation is 

X* = (1 ~ r 2 )NVar(y! . . . yjv) (15.2.13) 

where 

NVar(yi . . ,y N ) = ~ (15.2.14) 

- »=i 

For data with varying weights a it the above equations remain valid if the sums in 
equation (143.1) are weighted by l/<r?. 

The following function, fit, carries out exactly the operations that we have 
discussed. When the weights a are known in advance, the calculations exactly 
correspond to the formulas above. However, when weights a are unavailable, 
the routine assumes equal values of a for each point and assumes a good fit, as 
discussed in §15.1. . - 

The formulas (15.2.6) are susceptible to roundoff error. Accordingly, we 

rewrite them as follows: Define 
and 

Then, as you can verify by direct substitution. 



(15.2.17) 
(15.2-18) 



(15JL15) 
(15.2.16) 
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(15.2.19) 




(1512.20) 




(15.2.21) 




(15^22) 



'include <*&th.n> 
♦include ~nmtil.h" 

roid fitttlMt xD. float yO. int ndata, float sigQ, Int art, float *a, 
float *b, flomt *»iga. float •eigb, float *chi2. float *q) 

2^1* nLlal^T* *J^?i?'7 Cl - nda1:a3 indi **»' ****** deviations 
* h AV^?^1L^L ?\? straight line y = a + te by minimizing ^.Returned are 
^a^l^^^r^ ""^ntJes siga andsigo. the cW^quare ch±2. and the 
fr^T^f"^ P^^TNy q (that the fit would have x* this large or larger). If Bnrt=0 on 
mput. then the standard devotions are assumed to be unavailable: q is returned as 1 O and 
the Bormafization of ch±2 is to unit standard deviation on afl points? 

float gaaMqCfloat a, float x); 
lnt 1; 

float wt,t r 8Xoaa.BX-0.0.«y-0.0,rt2-0,0 r »» F «lgdait; * 
•b-O.O; 

** (m ^l i Accumulate sums ... 

for a»l;l<-ndata;l^) { _with weights 

a» wt; 
ax zCU*vt; 
^ sj +- yCtf*nt; 

> aloa { 

for a-l;lomi&ta;i++) < _or without weights, 

ax xfl]; 
ay*-yU); 



> 

axoea-ex/aa; 
If (anrt) < 

for a-l;l<-9data;l++) { 
t-Cx tl)~axoas)/aigti3 ; 
at2 t*t; 
^ »b t^Ci3/algCU; 

> ala* < 

for (l-l;l<-adata;i++) < 

t»x Cl3 -sxoea ; "* 
at2 *- t*t; 
*b t*yClJ; 

> * 

I? WAW Sofve foro.6, a., and o>. 

•a»C»y-axa<*W>/aa; * 

•Blgo-sqrt (1 . 0/st2) ; 



666 



Chapter 15. Modeling of Data 



♦chi2-0.0; Calculate X *. 

if . (mrt — 0) < 

for (i-l;i<«ndata;i++> 

*chi2 +- sqB(y[iJ-(*a>-(*b)*xC13>; 

•q«1.0; 

aigda*-Bqrt((*chi2)/(nilata.-2)) ; For unweighted data evaluate typt- . 

*oiga *- aigdat; «*S using ch±2, and adjust 

*8lgb *- eigdat ; the standard deviations. 

> else { 

for <l»i;i<-ndata;l++> 

*chi2 -H- SQR((yU3-(*a)-<*b)*xCi3)/«lgCl3>; 
*q=ganmq(0 . 6* <ndata-2> ,0 . 6* C*chi2> ) ; Equation (15.2.12). 



CITED REFERENCES AND FURTHER READING: 

Bevington. P.R. 1369. Data Reduction and Error Analysis for (he Physical Sciences (New Yortc 
McGraw-HB), Chapter 6. 



15.3 Straight-Line Data with Errors in Both 
Coordinates 

. If experimental data arc subject to measurement error not only in the jft's. bat also is 
the xi's, then the task of fitting a straight-line model 

y(x) = o + fex (153.1) 
is considerably harder. It is straightforward to write down the x* merit function for this case, 

*'<«■»> ~S w 

where <r xi and a v % are, respectively, the x and y standard deviations for the ith point Toe 
weighted sum of variances in the denominator of equation (153.2) can be understood both 
as the variance in the direction of the smallest x* between each data point and the Kne with 
slope 6, and also as the variance of the linear combination — a — oxt of two random 
variables Xi and y«, 

Varfo -a-bxi) = Varfo) + 6 3 \ar(«<) = <rj 4 + b*<rli = 1/w (1533) 

The sum of the square of N random variables, each normalized by its variance, is thus 
X* -distributed. ^ 

We want to ininimize equation (153.2) with respect to a and b. Unfortunately, the 
occurrence of & in the denominator of equation (153.2) makes the resulting equation for 
the slope &y? jdb = 0 nonlinear. However, the corresponding condition for the intercept, 
dx^/Oa = 0, is still linear and yields 

« » ^2 Wifa - 6x,)J ~ ( 1534 > 

where the uh's are defined by equation (1533). A reasonable strategy, now, is to use the 
machinery of Chapter 10 (eg., the routine brent) for niimmizmg a general one~dmiesisional 
function to mmimize with respect to 6, while using equation (153.4) at each stage to ensure 
that the minimum with respect to 6 is also minimized with respect to a. 



15.3 Straight-Line Data with Errors in Both Coordinates 667 




